RATE DETERMINATION FROM VECTOR OBSERVATIONS 


By 

Jerold L. Weiss 

Ball Space Systems Division, Boulder CO 

ABSTRACT 


N 9 3--ZA 

/•; yV '■> 7 


1 


Vector observations are a common class of attitude data provided by a wide variety of attitude sensors. Attitude 
determination from vector observations is a well-understood process and numerous algorithms such as the TRIAD algorithm 
[Shuster] exist. These algorithms require measurement of the line of site (LOS) vector to reference objects and knowledge 
of the LOS directions in some predetermined reference frame. Once attitude is determined, it is a simple matter to 
synthesize vehicle rate using some form of lead-lag filter, and then use it for vehicle stabilization. Many situations arise, 
however, in which rate knowledge is required but knowledge of the nominal LOS directions are not available. This paper 
presents two methods for determining spacecraft angular rates from vector observations without apriori knowledge of the 
vector directions. The first approach uses an extended Kalman filter with a spacecraft dynamic model and a kinematic 
model representing the motion of the observed LOS vectors. The second approach uses a "differential’' TRIAD algorithm to 
compute the incremental direction cosine matrix, from which vehicle rate is then derived. 


1.0 INTRODUCTION 

The development of small, light-weight, low-cost spacecraft that are reasonably complex requires an attitude 
determination and control system (ADCS) design in which the ADCS hardware supports more than a single operating mode 
or function. For example, a three axis magnetometer can be used in aiding reaction wheel desaturation and in support of 
coarse attitude determination; reaction control jets can be configured to provide both delta-V and control torques; and a 
horizon scanner that can be locked in place can support both spinning and 3-axis stabilized modes of operation. 

The primary motivation for the work reported in this paper came from applying this principle to the following 
problem : Given a zero-net momentum ADCS that includes a 4rt Sr sun sensor (used in a sun-seek safe mode), a 3 axis 
magnetometer (used in reaction wheel desaturation), and reaction control system (RCS) jets (used during orbit raising); 
provide a capability to detumble the spacecraft following large launch vehicle separation transients. Normally, this problem 
would be solved by adding a gyro package to the suite of ADCS hardware. The gyro data would be used to measure body 
rates that would be fed back to the RCS jets to detumble the spacecraft. However, to keep costs and weight down the gyro 
package is undesirable. The alternate solution that was developed was to use the existing hardware to measure the motion of 
the sun vector and the Earth's magnetic field vector and use this motion to (approximately) determine the vehicle's inertial 
rate. One way of accomplishing this would be to use an algorithm such as the TRIAD algorithm [Shuster] to compute the 
vehicle attitude and then determine rate with a lead-lag filter. However, this would require both a current vehicle ephemeris 
and an orbit propagator since the direction of the sun and magnetic field vector in a non-rotating earth-centered frame is 
required. The methods presented in this paper provide alternatives which do not require either an ephemeris or orbit 


propagator. 



A secondary motivation for this work comes from recent gyro-failure problems on the Hubble Space Telescope. 
Currently, the redundant gyros on HST are working and will probably be used to stabilize the spacecraft for shuttle-capture 
during the upcoming repair mission. If the redundant gyros were to fail, however, a method of stabilizing the spacecraft 
using only star tracker data may be needed. The work presented in this paper could be used to stabilize a spacecraft using 
star tracker data without performing the complex star identification process. 

The remainder of this paper is organized as follows. Section 2 develops a standard extended Kalman Filter (EKF) 
approach to the problem of estimating spacecraft rate using LOS vector observations without apriori knowledge of the LOS 
reference directions. Section 3 presents an alternative approach which is less accurate but computationally more efficient. 
This approach is based on the TRIAD algorithm and is called DTRIAD for "differential" TRIAD. Section 4 presents several 
simulation results using both methods and compares the rate determination accuracy and vehicle stabilization capability for 
the two methods. Finally, section 5 presents conclusions and suggestions for further work. 

2.0 EXTENDED KALMAN FILTER (EKF) 

The EKF approach applies the standard EKF equations [Gelb] to a system whose state is defined as the body- 
referenced angular rate vector, w, augmented with the body-referenced vectors pointing to the n reference objects that are 
potentially observable, e v i = 1 ... n. The next two sections develop the state and covariance propagation and update 
equations . For clarity of presentation in the following, assume that n = 2. 

2.1 State and Covariance Propagation 

The nonlinear state dynamics are derived from Euler's equations and the well known expressions for the motion of 
vectors in a rotating frame. They are : 

(!) " = 

(2) e, = -to x e i i = 1 ... n 

where N eit is the vector of external torques applied to the vehicle, / is the vehicle’s inertia matrix and h wU is the stored 
momentum in reaction wheels (/ includes the the inertias of the wheels; see [Wertz]). The implicit assumption in Eq. 2 is 
that the e, are fixed in a non-rotating inertial reference frame. 

At time the EKF propagation step requires that these equations be integrated from time to t^j using the 
updated estimate at time t * as a starting point. Numerous integration methods are available (Runge-Kutta, 
Forward/Backward Euler, Trapezoidal) each with its own merits (e.g., see [Hildebrand]). When the sample rate is fast and 
the magnitude of of to is small, the simplest method is Forward Euler or rectangular integration. This method will be used in 
Section 4.0 to demonstrate results. 

In addition to propagation of the state vector, the EKF must propagate the state covariance matrix. The typical 
approach is to assume small deviations of the state from the estimate between t k to t^+j. This implies that the linearized 
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equations of motion adequately represent the error dynamics and can therefore be used to form the Lyapunov equation that 
represents the dynamics of the covariance matrix. For this problem, the linearized dynamics are : 
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and where we have assumed a diagonal inertia matrix, I. The extension to nondiagonal inertia matrices and to n > 2 is 
straightforward. The 3x1 process noise vector, w, is assumed to be a zero mean white Gaussian noise process with 
covariance matrix Q. It is intended to represent the uncertainty in N ext but is also used to represent expected errors in the 
integration process for state propagation. 

Given the linearized perturbation equation (3), the covariance propagation equation is the standard Lyapunov 
equation : 


(4) 


dP/dt = AP + PA t + BQB t 


where P is the covariance matrix of perturbation state 5x and A and B are evaluated at x(k) . This can be solved by any 
number of methods. A common simple approach [Jazwinski] that works well in practice is to discretize (3) and use the 
discrete Lyapunov formulation, viz. 

(5) p-(k+l ) = F P + (k) F t + GQG t 


where F = I + AA, G = BA, P'(k+1 ) is the covariance matrix at time j before updating with measurements at time 
P + (k) is the covariance matrix at after updating with measurements at time //. and A = - ifc. Another simple approach 

is to integrate Eq. (4) using forward Euler integration. 
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2.2 State Update 

The second step in the EKF approach updates the state estimate based on measurements and requires an observation 
equation that relates the measurements to the states. The actual observation equation depends on the type of sensor being 
used. For example, a vector magnetometer directly measures the earth's magnetic field vector resulting in the simplest 
measurement equation 

(6) y = + v 

where v is the measurement error.* A two-axis sun sensor or a star tracker, however, produces two measurements that are, in 
general, a nonlinear function of e, (e.g., two components of the orthogonal projection of e, onto an image plane). In these 
cases, it is a straightforward matter to derive the general nonlinear observation equation 

(7) y = h(ei) + v 

and to determine the Jacobian matrix H = 8y/8e ( -. The standard EKF update equations are then 


(8) 

P*(k) = {p~(k)~ l +h t r~ 1 h] 

(9) 

K = P+(k) H t R 1 

(10) 

x(k) = x(k)+K[y-h(I i (k))] 


where T(k) is the result of integrating Eq. 's (1) and (2) forward from t k .j to t k starting at i(k-l), the e t (k) are the 
reference vector components of x(k), and R is the covariance matrix of v. The inverse formulation for the covariance 
matrix is used here to avoid numerical problems in the subsequent simulations. It is not the fastest implementation (e.g., see 
[Bierman]). Further work is needed to establish the best scaling and covariance update method for this problem. Note that 
the option of transforming the nonlinear measurements into a vector observation and then employing (6) as the measurement 
equation was not considered since the correct corresponding value of R is singular. When R is singular, P + can become 
singular, [Gelb] and K can not be calculated. If the transformation method is used, an approximate nonsingular R matrix 
should be used (e.g., /?,, = arcsine( angle error ) ). 


Strictly speaking, the earth s magnetic field violates the assumption that the vector £ j be fixed in a nonrotating inertial frame since the direction varies as the 
spacecraft orbits the earth. The resulting rate errors are insignificant in several applications. 
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2.3 Implementation Notes 

(1) The EKF as described above must be initialized with a state and covariance estimate. The method described in the 
next section is an excellent and simple way to obtain an initial state estimate. Further work is needed to derive the 
corresponding initial covariance estimate corresponding to this state-initialization method. 

(2) Choice of sample rate, A, depends on accuracy requirements and is influenced by the magnitude of the vehicle 
angular rate, co, and applied torques, N ext and by processing resources. 

(3) While Q is normally used to represent uncertainty in N eX f it should also account for the errors made in integrating 
the state equations (1) and (2). As with most EKFs, it is likely that some degree of tuning of the parameter Q will 
be needed. 

(4) The covariance update equation (8) as written requires inversion of a x 3n+3 matrix. As mentioned, this is 
not the most efficient propagation equation. The experience gained so far with simulation of this approach suggests 
that the covariance matrix P and hence the Kalman gain K converge to steady state values that are a function of the 
problem geometry (angle between the vectors e j and ^ 2 )* The rate of convergence appears to be dependent on 
vehicle rate. These results suggest that for constant problem geometry a constant gain EKF (which does not require 
on-line covariance propagation or update ) might provide satisfactory performance. 

(5) While it is clear that the angular rate vector is not completely observable when n = 1, the EKF equations can still be 
executed in this case. If a good estimate of co is initially available, then the result should be that the component of 
the rate estimate that is aligned with e^has a random walk error behavior. This is a convenient feature in cases 
where sensor dropouts occasionally occur. 

(6) Measurement biases and sensor misalignments will not affect rate determination accuracy (for this method or the 
one presented in the next section) since neither affects the change in LOS measurements from one time to the next. 

3.0 DIFFERENTIAL TRIAD APPROACH 

While the EKF approach described in Section 2 provides near-optimal solutions and can operate, with reduced 

accuracy, when n=l, the computational complexity of the algorithm can be relatively large for the current generation of 

spacecraft flight computers. The approach discussed in this section provides a simpler alternative. 

The idea for this approach comes from the basic strapdown equations of motion for the direction cosine matrix. A, 

[Wertz], 


( 11 ) 


A = £2(ft))A 


Assuming that a) is constant over the interval t £ to 


(12) 

A(t k+1 ) = e^ A A(t k ) 

(13) 

A' = A(t k+1 )A(t k ) T =e Q (“> A 
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which implies that knowledge of the direction cosine matrix representing the transformation from the vehicle attitude at time 
t k lo time t k+ ] (A") is sufficient for determining to. 

The TRIAD algorithm [Shuster, Wertz] is a well known method for determining the direction cosine matrix from 
vector observations assuming that the direction of the vector observations in some arbitrary reference frame is known. If this 
reference frame is the vehicle's attitude at time t k then it is clear that the TRIAD algorithm produces the matrix A' from 
which the vehicle body rate can be derived. Let e t m (k), i=l,2 be the measured vectors to two reference objects at time k. 
The "Differential TRIAD" (DTRIAD) approach to rate determination is as follows : 

(Step 1) Form the "reference" matrix A r from 

( 14 > = [>r (*) j (k)xe? ( k ) j er (k)x(e } m (k)xe? (*))] 

(Step 2) Form the "observation" matrix 

< 15 > = [«*(* + 1) ! e?(k + \)xe!?(k + \) | er(* + l)xUr(*+l)xe 2 m (* + l))] 

(Step 3) Compute A' = A m . 

(Step 4) Solve (13) for to. If the magnitude of Adi is sufficiently small , then the identity 
(16) A" -1 + £l(o))A 

makes this step straightforward. As the magnitude of A to increases, one can expand (13) in the standard power series and 
solve the polynomial equations that result; or one can reduce A until acceptable accuracy is achieved with (16). These steps 
arc repeated at every t k . Recognizing that A r (k) = A m (k-1) , we see that at each time t k , only steps 2 - 4 need to be 
repeated. 

The performance of the DTRIAD algorithm, like the performance of TRIAD, depends on problem geometry. The 
accuracy of the rate estimate in the direction which bisects the two reference vectors is the least accurate and can be 
estimated using standard methods. Like the EKF, performance will degrade with large values of Aco. Unlike the EKF, 
however, degradation of the estimate will also occur when large external torques are applied. Also, like the EKF, 
measurement biases and sensor misalignments will not affect the rate determination accuracy. 

4.0 SIMULATION RESULTS 

This section presents two simulation cases to demonstrate performance of the two algorithms. The first case 
assumes that the vehicle is tumbling at a relatively high rate (say -10 deg/sec resulting from an open loop despin operation 


280 



following an SRM firing) and that continuous measurement of the Earth's magnetic field and the solar vector are available. 
For this case, we show both the accuracy of the two rate determination methods and the resulting tumble recovery capability 
assuming a feedback control law using reaction control jets. The second case assumes that the vehicle is slowly tumbling 
(say ~ldeg/sec as the result momentum accumulation due to external disturbances and wheel run-down following the failure 
of primary attitude control gyros) and that continuous tracking of stellar objects from a star tracker are available. In this 
case, we show the rate stabilization capability using both methods and reaction wheels for control. In all cases the EKF is 
initialized by adding a small error to the true state. In practice, the EKF is initialized using the DTRIAD algorithm. All the 
simulations presented herein are full nonlinear 3 degree of freedom simulations that use 6’th-order Runge-Kutta integration 
to solve Eq. 1 along with the standard quaternion kinematic equations. 

4.1 Hiph-rate Tumble Recovery using S un Sensor and Magnetometer Data 

For the purposes of this simulation, we assume that both the sun sensor and the magnetometer provide 
measurements of reference vectors (i.e., Eq. (6) applies). It is assumed that the sun vector measurement is corrupted by a 
constant bias corresponding to a 1 degree rotation of the true sun vector and a white noise whose magnitude corresponds to 
0.25 degrees error (i.e., R = (.25*tc/ 180) 2 *E, where E is the identity matrix). The magnetic field vector measurement is also 
corrupted by a 1 deg rotation bias and a white noise with a standard deviation of 9 nano Tesla (about 1 part in 25 for low 
earth orbits). The sensor errors correspond to commonplace sensor capabilities. The 1 degree bias error in the 
magnetometer vector includes errors due to electronic bias and misalignments of the sensor. The sun sensor bias error 
includes electronic bias, misalignments and Earth albedo effects (albedo effects are significant only for non-imaging cosine- 
law sun detectors ). 

Figure 1 shows the body rate profile and EKF-estimated rates for a vehicle with diagonal inertia matrix I = 
diag([156 85 186]) kg-m 2 and a 10Hz update rate (A=0.1). Figure 2 shows the corresponding EKF and DTRIAD rate errors. 
The EKF produces superior results due to the inclusion of a dynamic model which allows filtering to take place. 



Fig. 1 - Actual and EKF-Estimated Body Rates; 
N ext = 0, A=0.1 


* This result is somewhat artificial since a dynamic model can be used to form a simple filler (Kalman or otherwise) with the DTRIAD rate estimates used as 
"measurements". 
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Figure 3 shows the result when the EKF-determined rates are fed back to a reaction control system to stabilize the 
vehicle. In this simulation, it is assumed that the control system operations at a sample rate of 2Hz and that the reaction 
control jets are pulse-width-modulatcd using the standard impulse-matching algorithms. The thruster torques available in 
each axis are 1.1 Nm in roll 0.86Nm in pitch and 0.23 Nm in yaw. The jet pulse times are quantized in 20msec increments. 
Figure 4 shows the feedback response using the DTRIAD algorithm with all other parameters the same as for Fig. 3. The 
EKF is substanitally better than DTRIAD in this case because it includes a dynamics model that allows it to predict 


spacecraft motion based on applied RCS jet torques. 




Fig. 3 - Actual and Estimated Rates during 
Tumble Recover EKF Feedback. A=0.5, 

N ext = 


Fig. 4 - Actual and Estimated Rates during 
Tumble Recovery with DTRIAD Feedback. 

A =0.5, N ext = g(0. 






4.2 I .nw Rate Sta bilization using Star-Tracker Data 

In this simulation we assume that the star trackers provide reference vector measurements (Eq. 6 applies again). 
Two "guide star" measurements are used. They are separated by 5 deg and corrupted by a 6 arcsec noise-equivalent angle 
(typical values). For these simulations, the spacecraft inertia properties correspond to a large platform (I = 
diag([1000;1000;2000]) kg-m^) and control is provided by large reaction wheels (0.8Nm torque capability). In all 
simulations a 2 Hz sample rate is used and the initial total angular rate is 0.5 deg/sec. 

Figure 5 shows the magnitude of the vehicle rate and the rate errors using the DTRIAD algorithm. The 0.5 deg/sec 
rale error is damped in about 50 seconds and the RMS residual rate is between 0.0001 and 0.002 deg/sec per axis. Figure 6 
shows the rate magnitude and rate errors using the EKF. The damping time constant is the same as DTRIAD (the control 
gains are the same) and the RMS residual rate is between 0.00003 and .0002 deg/sec, much smaller than DTRIAD. The 
smaller residual rate is due to the better rate prediction accuracy which is due to the incorporation of applied torques in the 
EKF rate prediction algorithm. 
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Fig. 5 Actual Rate and DTRIAD Errors for Star Tracker Simulation; 

A = 0.5 and N ext = geo. 
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Fig. 6 - Actual Rate and EKF Errors Errors For Star Tracker Simulation; 

A = 0.5 andN ex( = geo. 


5.0 CONCLUSIONS 

Two methods of computing the body-referenced angular rate vector from vector observations without knowledge of 
the LOS vector directions in inertial space have been presented. The EKF approach is based on spacecraft dynamics models 
and the rotational kinematics of the observed reference vectors. The DTRIAD approach uses the TRIAD algorithm to 
compute incremental direction cosine matrices from which approximate rate is derived. The EKF is computationally more 
complex but produces superior estimates, especially in the presence of large control torques. Further work is needed to 
provide analytic bounds on accuracy for both methods, to determine the best covariance update implementation, and to 
determine the feasibility of using a constant gain EKF. 
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